** Installation of windows
window menu append submenu "stUser" "&Cluster RCTs"
window menu append item "Cluster RCTs" "Sample size for cluster RCTs (&clustersampsi)" "db clustersampsi"
window menu refresh

egen cluster1 = group(agency) if sample==1

mat B = J(`: word count $all_DV',3,.)
mat rownames B = $all_DV
mat colnames B = "MDE" "Beta" "SE"

xtreg business_w5 if sample==1 & time==0, i(cluster1) re
loc m  =round(e(g_avg))
loc k  =round(e(N_g)/2)
loc icc=e(rho)

loc A = 1+(((`m'/2)-1)*`icc')
loc r = ((`m'/2)*`icc')/`A'

sum business_w5 if sample==1 & treatment==0 & time==0
loc con_m = r(mean)
loc con_sd = r(sd)
 
** Code by Hemming and Masch considering baseline correlation for a continuous variable 
clustersampsi, detectabledifference mu1(`con_m') sd1(`con_sd') rho(`icc') alpha(0.05) beta(0.8) m(`m') k(`k') base_correl(`r')
loc rdd = r(DD)

mat B[1,1] = round(`rdd', .001)

qui reg business_w5 i.treatment $ControlsB if time==1, vce(cluster agency)

mat B[1,2] = round(_b[1.treatment], .001)
mat B[1,3] = round(_se[1.treatment], .001)

local i 1
foreach v of global all_binary {

	local ++i
	xtreg `v' if sample==1 & time==0, i(cluster1) re
	loc m  =round(e(g_avg))
	loc k  =round(e(N_g)/2)
	loc icc=e(rho)

	loc A = 1+(((`m'/2)-1)*`icc')
	loc r = ((`m'/2)*`icc')/`A'

	su `v' if sample==1 & treatment==0  & time==0
	loc con_m = r(mean)
	
	** Code by Hemming and Masch considering baseline correlation for a continuous variable 
	clustersampsi, binomial detectabledifference p1(`con_m')  rho(`icc') alpha(0.05) beta(0.8) m(`m') k(`k') base_correl(`r')
	
	loc rdd = r(DD_pos)

	mat B[`i',1] = round(`rdd', .001)

	reg `v' i.treatment $ControlsB if time==1, vce(cluster agency)
	
	mat B[`i',2] = round(_b[1.treatment], .001)
	mat B[`i',3] = round(_se[1.treatment], .001)
}

preserve
	svmat B, names(col)
	keep if Beta!=. 
	keep MDE Beta* SE*
	gen 	id=_n
	
	gen seUp90 = Beta + 1.645*SE
	gen seLow90 = Beta - 1.645*SE

	gen seUp95 = Beta + 1.96*SE
	gen seLow95 = Beta - 1.96*SE

	gen seUp99 = Beta + 2.576*SE
	gen seLow99 = Beta - 2.576*SE

	gen MDE_neg = MDE*(-1)
	
	
	la def variable 1 "Bus. days" 2 "In time" 3 "Complete" 4 "Open" 5 "Extension" 6 "0-15"  7 "0-25" 8 "15" 9 "25" /*10 "14/15/16" 11 "24/25/26" 12 "0-13" 13 "17-23"*/
	la val id variable
	
	tw (bar MDE id if id>=2, color(midgreen%30))(bar MDE_neg id if id>=2, color(midgreen%30))(rcap seUp99 seLow99 id if id>=2, lcolor(ltblue%70))(rcap seUp95 seLow95 id if id>=2, lcolor(teal%50)) (rcap seUp90 seLow90 id if id>=2, lcolor(navy))(scatter Beta id if id>=2, color(cranberry)), xtitle(" ", size(small)) xlabel(2 3 4 5 6 7 8 9 /*10 11 12 13*/, valuelabel labsize(medium)) ylabel( -0.24 "-24" -0.20 "-20" -0.16 "-16" -0.12 "-12" -0.08 "-8" -0.04 "4" 0 "0" 0.04 "4" 0.08 "8" 0.12 "12" 0.16 "16" 0.20 "20" 0.24 "24", angle(horizontal) labsize(medium)) ytitle("Percentage points (pp)", size(medlarge)) yline(0, lp(dash)) graphregion(fcolor(white) lcolor(white) ifcolor(white) ilcolor(white) margin(right)) legend(order(1 "MDE" 6 "ATE") pos(3) r(2) size(small)) name(pp, replace) 
	
	tw (bar MDE id if id==1, color(midgreen%30))(bar MDE_neg id if id==1, color(midgreen%30))(rcap seUp99 seLow99 id if id==1, lcolor(ltblue%70))(rcap seUp95 seLow95 id if id==1, lcolor(teal%50)) (rcap seUp90 seLow90 id if id==1, lcolor(navy))(scatter Beta id if id==1, color(cranberry)), xtitle(" ", size(small)) xlabel(1, valuelabel labsize(medium)) ylabel(-5(1)5, angle(horizontal) labsize(medium)) ytitle("Number of days", size(medlarge)) yline(0, lp(dash)) graphregion(fcolor(white) lcolor(white) ifcolor(white) ilcolor(white) margin(right)) legend(off) name(bd, replace) fxsize(25) fysize(100)
	
graph combine bd pp, holes(3) rows(2) imargin(0 2 0 0) graphregion(fcolor(white) lcolor(white) ifcolor(white) ilcolor(white) margin(right)) name(coeff, replace)
graph display coeff, xsize(20) ysize(12) margin(left) 

graph export "$Figures/Figure_C1.pdf", as(pdf) replace

